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In numerical simulations of black hole binaries, Pretorius and Khurana [Class. Quant. Grav. 24, 
S83 (2007)] have observed critical behaviour at the threshold between scattering and immediate 
merger. The number of orbits scales as n ~ — 7 In \p — p* | along any one-parameter family of initial 
data such that the threshold is at p — p„. Hence they conjecture that in ultrarelavistic collisions 
almost all the kinetic energy can be converted into gravitational waves if the impact parameter is 
fine-tuned to the threshold. As a toy model for the binary, they consider the geodesic motion of 
a test particle in a Kerr black hole spacetime, where the unstable circular geodesies play the role 
of critical solutions, and calculate the critical exponent 7. Here, we incorporate radiation reaction 
into this model using the self-force approximation. The critical solution now evolves adiabatically 
along a sequence of unstable circular geodesic orbits under the effect of the self-force. We confirm 
that almost all the initial energy and angular momentum are radiated on the critical solution. Our 
calculation suggests that, even for infinite initial energy, this happens over a finite number of orbits 
given by n m ~ OAl/n, where r\ is the (small) mass ratio. We derive expressions for the time spent 
on the critical solution, number of orbits and radiated energy as functions of the initial energy and 
impact parameter. 
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I. INTRODUCTION 

Motivated by speculation about the formation of black 
holes in high-energy particle collisions, there has recently 
been interest in binary black hole mergers with large ini- 
tial boosts; see [l[ for a review. A critical surface in 
phase space separates initial data which lead to merger 
from data which scatter at first approach, thus defining a 
threshold of immediate merger in the space of initial data. 
(Note that "scatter at first approach" here includes data 
which merge at a subsequent approach.) 

Pretorius and Khurana 0] have numerically evolved 
initial data near this threshold for equal mass, non- 
spinning black hole binaries. They find that generic 
smooth one-parameter families of initial data intersect 
the threshold once, consistent with it being a smooth hy- 
persurface in phase space. Moreover, they find that the 
number n of orbits before either merger or flying apart 
increases with fine-tuning to the threshold, with n scaling 
approximately as 



7 In \p — p* I + constant, 



(D 



where p is any smooth parameter along the family of ini- 
tial data, and its value at the threshold of immediate 
merger. The additive constant depends on the family 
and on how p is defined, but the dimcnsionless critical 
exponent 7 is independent of the family (for the reasons 
reviewed in Sec. [TTJ) . In Q, where initial data with rela- 
tive boosts of around 0.22 are evolved, ((T|) is found to be 
a good fit for 1 < n < 4 with 7 = 0.35 ± 0.03. 

In this regime, a fraction of 1 to 1.5% of the rest mass 
is radiated per orbit 0. While their initial data are rest 
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mass-dominated, Prctorius and Khurana speculate that 
in large-boost initial data fine-tuned to the critical im- 
pact parameter, most of the initial kinetic energy, and 
hence most of the total energy, can be turned into grav- 
itational radiation over a large number of orbits. This is 
quite different from the situation for either large Q or 
small Q impact parameter. Hence in a particle physics 
context, there is a small cross section (corresponding to 
fine-tuned impact parameter) in which the gravitational 
wave energy emitted is much larger than in a generic 
collision. 

In a subsequent paper usinga different numerical code, 
Sperhake and collaborators [5( evolved initial data with 
boosts around 1.5 and 2.9. The critical exponent 7 is 
estimated as 0.2, and the scaling law is observed for the 
range 0.5 < n < 1.5. No data for critical scaling of the 
energy are given but in these more relativistic collisions, 
the fraction of total energy radiated is reported to be 
as high as 35% for the higher boost. Comparison with 
the results of Q suggests that 7 and the energy radiated 
depend on the initial boost (initial energy). 

Pretorius and Khurana @] note that what they observe 
in full numerical relativity is similar to "zoom-whirl" be- 
haviour Q in the timelike geodesies of test particles in a 
black-hole spacetime when the impact parameter is fine- 
tuned to its critical value, now characterising the thresh- 
old of immediate capture. Hence they propose equatorial 
orbits of a test particle on a Kerr spacetime (with mass 
and angular momentum set to the total mass and angular 
momentum of the binary) as a toy model for the equal 
mass binary. They point out that the critical solutions 
mediating this behaviour are the unstable circular orbits, 
and use this to calculate the critical exponent for this toy 
model. We obtain their results as a limiting case of ours 
in Sec. |VC] below. 

The presence of critical phenomena suggests the ex- 
istence of a critical solution defined by having a sin- 
gle growing perturbation mode and usually also charac- 
terised by a symmetry. Because energy is lost through 
gravitational radiation, the critical solution cannot be 
stationary, and because the black holes have rest mass it 
cannot be self-similar. We suggest that the critical so- 
lution is an unstable "as circular as possible" orbit that 
in the limit where one of the two objects becomes a test 
particle is indeed one of the unstable circular orbits con- 
sidered by Pretorius and Khurana. 

In this paper, we extend their analysis to binaries with 
large but finite mass ratio, and include radiation reac- 
tion in the self-force approximation. We restrict to the 
case where the large black hole is Schwarzschild, but our 
methods are in principle applicable to Kerr. This math- 
ematical approach allows us to clarify the nature of the 
critical solution as adiabatically stationary and, based on 
this, compute the orbit and radiation as a function of the 
initial data. 

In the large mass ratio regime, we confirm the con- 
jecture of Pretorius and Khurana that almost all the ki- 
netic energy of the binary (which in turn is almost all the 



total energy for ultrarelativistic collisions) is radiated if 
the fine-tuning of the initial data is sufficiently good. We 
find, however, that the number of orbits remains finite 
as fine-tuning is improved, while ((T|) holds only for suffi- 
ciently small 77 and \p — | . 

We begin by reviewing the general mathematical ideas 
behind critical phenomena in dynamical systems lan- 
guage in Scc.|Hj For reference, and to establish notation, 
we review timelike geodesies in Schwarzschild spacetime 
in Sec. HID We construct the critical solution in the self- 
force approximation in Sec. IIV1 and in Sec. |V] we inves- 
tigate its perturbations and from these we derive expres- 
sions for the the total number of orbits and total energy 
radiated as functions of the initial energy and impact pa- 
rameter. Sec. I VII summarises our results and gives an 
outlook on the comparable mass case. 

Throughout the paper ~ denotes equality up to sub- 
leading terms, and ~ equality up to sublcading terms and 
an overall constant factor, and := denotes a definition. 



II. DYNAMICAL SYSTEMS IDEAS 

We briefly review type I critical phenomena and the 
calculation of the related critical exponent, in the lan- 
guage of abstract dynamical systems. Let <f> represent a 
point in phase space. This could be either a vector of 
variables or a vector of fields at one moment of time. Let 
4>{t) be a trajectory in phase space that represents a so- 
lution. For the basic critical phenomena picture, we do 
not need to distinguish between field theories and finite- 
dimensional dynamical systems (and so do not write any 
x-dependence). This will later allow us to approximate 
a field theory (general relativity) by a finite-dimensional 
dissipative dynamical system (particle orbits with self- 
force) . 

Assume that at late times there are two qualitatively 
distinct outcomes, such as scattering and plunge in our 
system. There must then be a hypersurface in phase 
space separating these two basins of attraction, called the 
critical surface. As solution curves in phase space cannot 
intersect, the critical surface must be a dynamical system 
in its own right. Assume that the critical surface, consid- 
ered as a dynamical system, has an attractor, called the 
critical solution 0*. As 0* is a fixed point (independent 
of t), its linear perturbations must be a sum of modes of 
the form e Xit <fii, with <f>i also independent of t. As <f>„ is 
an attractor within the critical hypersurface, which itself 
is a repeller, must have precisely one growing mode, 
pointing out of the critical surface, for some (real) Ao > 

h 

Consider a one-parameter family of initial data, with 
parameter p, such that this family intersects the critical 
surface at p = p* . The evolution of the precisely critical 
initial data with must lie in the critical surface and 
hence will find Consider now initial data near the 
critical surface, with p — p* sufficiently small so that <f>(t) 
passes close to 0* during the evolution, and we can use 
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perturbation theory about 4>* ■ There is then a range of t 
where the decaying modes can be already neglected and 
the growing mode 4>q is still small enough for perturba- 
tion theory to hold, leading to the approximation 



<j>{t) ~0*+Ci(p-p*)e A °V, 



(2) 



where C\ is some constant that depends on the particular 
one-parameter family. Now define by 



(3) 



where C'2 is a fairly arbitrary constant, representing a 
reference deviation from the critical solution where it be- 
comes apparent on which side of the critical surface the 
solution is going to end up. One then finds 



t.^-rin|p-p,|+C 3 



(4) 



where C3 depends on C\ and C2, and T = 1/Ao is called 
the critical exponent. 

Essentially the same picture holds if (f>* is not a fixed 
point but a limit cycle (periodic in t). Then 4>q is also 
periodic in t and so is C3, and this leads to a modulation 
of the scaling law periodic in In \p — p* | . This generalisa- 
tion is not relevant for our application, but might be for 
extreme-mass ratio orbits in Kerr not in the equatorial 
plane, or more generally for binaries with spin. 

However, another generalisation is relevant here. As 
we shall see, in the geodesic toy model of Pretorius and 
Khurana the critical solution is not an isolated attractor, 
but instead there is a line of critical solutions (each of 
which represents an unstable circular orbit, and has pre- 
cisely one growing mode) . In the self- force approximation 
to extreme mass ratio binaries, by contrast, these merge 
into a single critical solution, which can be approximated 
as an adiabatic motion along the line of unstable circu- 
lar geodesic orbits, and so is slowly time-dependent. The 
general dynamical system picture above remains again 
basically unchanged, except that in the first case </>*, 4>o 
and Ao form a one-parameter family (parameterised for 
example by the energy of the orbit), and in the second 
case they become slowly varying functions of t. 



III. TEST PARTICLES ON SCHWARZSCHILD 
SPACETIME 

A. Equations of motion 

To establish notation, we review the trajectories of test 
particles, modelled as geodesies, on Schwarzschild space- 
time. Let t, r, 0, if be the usual Schwarzschild coordi- 
nates, so that the Schwarzschild metric is 



ds 2 



l--\df 

r 



2 

1 - - 

r 



dr 2 



+r 2 {d0 2 + sin 2 ddip 2 ), 



(5) 



and let a dot denote the derivative with respect to the 
proper time r of the particle. Throughout this paper we 
use units such that Newton's constant G, the speed of 
light c and the mass M of the Schwarzschild spacetime 
are all unity. Let u a = x a = dx a jdr be the 4- velocity of 
the particle. Without loss of generality let the orbit be 
in the plane 9 = ir/2. Define 



E := 
L := 



dt 



U a , 



(6) 
(7) 



representing the energy and angular momentum per rest 
mass of the test particle. Retaining the convention c = 
G = 1, but restoring M, and with m the mass of the test 
particle, its physical energy and angular momentum are 
mE and mML. As d/dt and d/d(p are Killing vectors 
of the Schwarzschild spacetime, E and L are conserved 
quantities in the sense that 



E = 0, 
L = 

along geodesic orbits. Hence we obtain 



t= 1- 



E. 



<i>- L 



(8) 
(9) 



(10) 
(11) 



By definition, E > for any timelike geodesic for r > 2, 
where d/dt is future timelike. Without loss of generality 
we shall also assume L > 0. The normalisation condition 
u a u a = — 1 can then be expressed as 



where 



E 2 



V{L,r) 



r 2 + V(L,r), 



1 



1 



Geodesies on Schwarzschild obey 

u b V b u a = x a + T a bc x b x c = 0, 



(12) 



(13) 



(14) 



but (|8I9[) together with the normalisation condition ([12")) 
can be used to reduce this 6-dimcnsional dynamical sys- 
tem to a 3-dimensional one in the variables (r, f, L), with 
the equations of motion 



r = -\vAL,r), 
L = 0, 



(15) 
(16) 



and where E(r,f,L) defined by (TT2")) is an integral of the 
motion. (Here and in the following, commas denote par- 
tial derivatives.) Alternatively, as long as r does not 
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FIG. 1. For Sec. IIII At Illustration of the effective po- 
tential (|13[) for different values of L. From above, the four 
full lines illustrate the cases L > 4, VT2 < L < 4, L — \/l2 
and L < yl2. Note that all these curves approach V — V 1 
as r — > oo. - For Sec. IIV At The dashed line repre- 
sents the critical solution. Its almost vertical segment rep- 
resents the adiabatic phase. As illustrated by the beads, it is 
given by r(t) = r_(L(t)), E 2 (t) = V- (£(*)), with L(t) deter- 
mined by radiation. The horizontal segment represents the 
plunge phase. The two dotted lines represent orbits com- 
ing in from infinity at the critical impact parameter for dif- 
ferent E, and asymptoting to the critical solution. Radia- 
tion reaction has been neglected in the approach and plunge 
phases, while r has been approximated as zero in the adia- 
batic phase. With radiation reaction, to order r\ the approach 
and plunge phase curves bend slightly downwards, the adia- 
batic curve is slightly further to the right, and the kinks are 
slightly smoothed out. [The beads are the maxima of the po- 
tential for L — 5, 4.6, 3.9, \/l2, and the lowest curve (without 
a maximum) corresponds to L — 3.] 
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FIG. 2. For Sec. IIII Bt The space of geodesies on 
Schwarzschild, labelled by angular momentum per unit mass 
L and energy per unit mass (squared) E 2 of the particle. The 
bottom right shaded region contains bound orbits that os- 
cillate between a minimum and maximum radius. Its lower 
edge is given by stable circular orbits. The left shaded region 
contains orbits that come in from infinity with low impact 
parameter and plunge. The top right shaded region contains 
orbits that come in from infinity with high impact parameter 
and scatter. The line dividing these two regions is therefore 
the threshold of immediate capture, given by E 2 = V-(L) 
(which characterises unstable circular orbits). The white re- 
gion corresponds to orbits that emerge from the white hole 
horizon and plunge back into the black hole, and which are 
not of interest to us here. — For Sec. IIV At In the self-force 
approximation, the critical solution moves adiabatically along 
the line of unstable circular geodesic orbits, as indicated by 
the dashed line and arrow. The beads and arrow are as in 
Fig- El (The approach and plunge phases seen in Fig. [1] are 
essentially in the r direction, which is suppressed here.) 



change sign, the system can be written in the variables 
(r, L, E) in the form 



f = ±v/£ 2 - V(L,r), 
L = 0, 
E = 0, 



(17) 
(18) 
(19) 



which now shows both integrals of the motion explicitly. 
In either case, the variables (t, <p) are evolved by (|1Q|11|) 
and play no dynamical role, while the remaining compo- 
nent of (|14[) is consistent with our assumption 9 = 0. 

The shape of the effective potential V(L, r) for different 
ranges of L is illustrated in Fig. [1] (The dashed and 
dotted curves should be ignored for now.) As V(L, 00) = 
1 for any L, orbits with E > 1 are unbound. For L > 
\/l2, the effective potential has a maximum at r = r_ 



and a minimum at r = r + , given by 

L(L±L 12 ) 



where it takes values 



V±:=V(r±) = 



2(2L±L 12 ) 2 
9L{L±L 12 ) 

Here we have introduced the shorthand 

L u := ^/L 2 - 12. 



(20) 



(21) 



(22) 



The minimum and maximum give rise to a stable and 
an unstable circular orbit. These merge for L = ^/\2, at 
r + = r_ = 6. 

For trajectories coming in from infinity, and hence with 
E > 1, a natural parameter of the initial data is the 
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impact parameter b > 0, defined by 

b := lira r sin \p(f) — <^(oo)|. 

r— >oo 

With 



Perturbations of the critical solution 



we find 



dr r 2 ^jE 2 - V(L,-, 



2 L 2 
E - 1= P- 



(23) 
(24) 

(25) 



We define the (instantaneous) orbital angular fre- 
quency 

o dtp <p [ 2\ L V L 

For both the stable and unstable circular geodesic orbits, 
we have 

E = , r ~ 2 ., (27) 



L = 



7^~3)' 
r 



(28) 
(29) 

We will be particularly interested in the ultrarclativistic 
limit £> 1, where 



L~ 3(r-3)~ 1/2 ~ 3V3E, 



B. Critical solution 



(30) 

(31) 
(32) 



For simplicity of presentation, in the following we dis- 
cuss only geodesic orbits with E > 1 and which initially 
have r < and r > 7'_(L), coming in from large r. If 
the impact parameter is fine-tuned so that E 2 = V~(L), 
the orbit asymptotically approaches the unstable circu- 
lar orbit r = r_(L), without ever plunging or moving 
back out. Orbits with the same L but smaller E scatter, 
while orbits with the same L but larger E plunge. Hence 
the critical surface in the space of initial data parame- 
terised by (r, E, L) with r < is formed by the surface 
E 2 = V-(L), r > r_(L). The space of geodesic orbits is 
illustrated in Fig. [5] 

We define the critical impact parameter b* (E) by E 2 = 
V-(L = L(E, 6*)) and find 



K(E) = 



V / 27E' 1 - 36E 2 



E(9E 2 - 8) 3 / 2 



(33) 



V2(E 2 - 1) 
which has the well-known limits 

fe*(oo) = 3-\/3 

and = oo (every particle released from rest at in- 

finity falls into the black hole). 



(34) 



The linearisation of the dynamical system (|15I16[) 
about the critical solution r = r_ (L) is 

5r = ~V >rr -8r - \v rL -5L, (35) 

SL = 0, (36) 

where V rr - and V^ r L- are V rr and V^ r L evaluated on 
the critical solution r = r_(L). We can decouple these 
equations as 



X = 1 X, 

SL = 0, 



where 



Sr 



Vrl 



V, 



r := ( --v rr . 



-SL, 



-1/2 



LV 2 {L~L 12 ) 



1/2 
12 



(37) 
(38) 

(39) 
(40) 



With E 2 = V-(L) on the critical solution, we can con- 
sider r also as function of E, but we have not been able 
to express T(E) in closed form. 

From (|35I36[) we see that the critical solution has ex- 
ponentially growing and decaying modes x — exp(±r/r) 
with SL = 0, and a neutral mode SL = 1 with x = 0. 
Hence it has precisely one growing mode, as required of 
a critical solution. However, in the geodesic approxima- 
tion we do not have a critical solution which is global 
attractor within the critical surface, but rather a line of 
critical solutions, to which the neutral mode is tangent. 



IV. CRITICAL SOLUTION IN THE 
SELF-FORCE APPROXIMATION 

A. Qualitative discussion 

Our aim is to replace the model of the binary merger 
as a test particle on Schwarzschild spacctime with a more 
realistic model that in particular incorporates radiation 
reaction, but which still approximates the Einstein equa- 
tions as a finite-dimensional dynamical system by inte- 
grating out the radiation. 

One such model is the extreme mass-ratio approxima- 
tion, where the loss of energy and angular momentum are 
calculated to leading order in 77 := m/M, where M = 1 is 
the mass of a large black hole and m that of the smaller 
object, modelled as a point particle. We shall use an adi- 
abatic approximation, where the particle is considered to 
move on a geodesic with E and L now slowly decaying, 
and a force term appearing in consequence on the right- 
hand side of the geodesic equation. We can then still use 
the effective potential picture to write the equations of 
motion, where E(t) and L(t) now evolve under the effect 
of the radiation reaction. 
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In short, the critical solution sits approximately at the 
maximum of the potential, with E{t) 2 = V-{L(t))+0(rj) 
and r(t) = r_(£(i)) + 0{rj), while that maximum it- 
self evolves under radiation reaction, starting with E ~ 
L 2 /27 and r ~ 3 as L — > oo and ending with L ~ s/V2, 
E ~ 8/9 and r ~ 6 at plunge. Our self- force calculation 
suggests that infinite E and L occur at finite r, i and ip in 
the past, but as we discuss below, we cannot consistently 
reach E = oo within the self-force approximation. 

Any initial data that are perfectly fine-tuned to the 
threshold of immediate merger, with arbitrarily large E, 
join the critical solution, which acts as an attractor in 
the critical surface. This is indicated schematically in 
Figs. Q] and [H where the dotted lines with arrows in- 
dicate two fine-tuned solutions approaching the critical 
solution (dashed) and remaining on it until and through 
the plunge. Note that the line of critical fixed points in 
the test particle picture has become a single, adiabati- 
cally evolving critical solution. 



B. Equations of motion 

In the extreme mass ratio limit, the larger object in the 
binary is approximated as a Schwarzschild black hole, 
and the lighter one as a particle moving in this back- 
ground spacctime, with the equation of motion 



b V b u a = F a , 



(41) 



where u & V&u a = describes a geodesic on the back- 
ground Schwarzschild spacetime, and F a is the self-force 
per particle rest mass, which is proportional to the mass 
ratio r\ to leading order. (We use terminology customary 
in the self-force literature where it is really mE, mML 
and mF a that have dimensions of energy, angular mo- 
mentum and force.) V a is the covariant derivative in the 
background spacetime, and indices are moved implicitly 
with the background metric. 

As before, we define E and L by (fTU)) and pT|) . and we 
impose the normalisation u a u a = — 1 in the form (|12p . 
The time derivative of u a u a = — 1 gives the usual con- 
straint on any 4-force, namely 



u a F a = 0. 



(42) 



From (|41j) . we now have the 3-dimensional dynamical 
system 



1 



V >r (L,r)+F r 



L = F ia 



(43) 
(44) 



E can be read off as E(r,r,L) from (fl"2"j) . or can be 
evolved using the third component of (|4"Tj) , which can 
be written as 



These are compatible because of (|4"2"j) . which can be writ- 
ten as 



rF r + E(F t + SIF V ) = 0, 



(46) 



where is defined by (ffijf . ip and t are evolved via (fTTj) 
and ([TOll. 



C. Adiabatic expansion 

By definition, the critical solution is the one that is 
balanced between plunging and running out to infinity. 
In the absence of radiation reaction, this clearly means 
it is a circular geodesic sitting on top of the potential. 
In the presence of radiation reaction, this definition be- 
comes teleological: the critical solution hesitates as long 
as possible between plunging and running out to infinity. 
In practice we use adiabaticity to define the critical so- 
lution as being as circular or as stationary as possible: r 
should vary over a radiation reaction timescale, that is, 
r should be proportional to the mass ratio r] to leading 
order in the self-force expansion. 

We formalize the adiabatic expansion through the se- 
ries ansatz (slow-time expansion) 



k=i 

oo 

M r ) =xy>fc(*)> 

fc=0 
oo 

oo 



k=0 



where we have defined the slow time 



T)T, 



(47) 
(48) 
(49) 
(50) 

(51) 



the suffix * denotes the critical solution, and the normal- 
isation condition (|12[) is obeyed order by order by virtue 
of (|4l))) . To leading order, 



hO(r ? 3 ), 



(52) 
(53) 
(54) 
(55) 



where / := df/dr as before, and /' := df/dr. Substitut- 
ing this into (|43l44p . we obtain to leading order [O(l) in 
O, and 0(77) in 0] 



E 



-Ft. 



(45) 



= V : r(L o ,r o ) 
L' = F 



r Q = r-(L ), 



¥>1 ' 



(56) 
(57) 
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and this allows us to integrate r and L . Hence, 
Fiij) = Fi( r o(T)) is the first-order self-force on a cir- 
cular geodesic orbit of radius ro, which is known Q. To 
next order, 

= -i[^ rr _(i ,r )ri + V; Pi (L ,ro)ii] + *T,(58) 

i'j = F v2 (59) 

formally give us ri and L\, but we cannot use this to 
calculate the critical solution beyond (ro,io). ^ i s 
the second-order self-force. It comprises terms that are 
quadratic in metric perturbations and corrections due to 
the fact that the critical solution is not exactly geodesic. 
Both have not yet been calculated. Note, however, that 
([551) can be rewritten as 

= -\v r (L Q + ijii.ro + vri) + VF[ + 0(r, 2 ), (60) 

and so its physical meaning is that the critical solution is 
slightly off the peak of the geodesic potential, held there 
(to this order) by the conservative part F[ of the self- 
force. 
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FIG. 3. Energy loss E/rj as a function of z := 1 — 3/r, for 
a circular geodesic orbit in the self-force approximation, and 
two closed-form fits. The fitting ansatz is given by Eq. (J6TJ) 
with a either fitted (solid line), or fixed at a — 3/2 (dashed 
line). The self-force data (points) are given in Appendix |A"1 
The insets show the relative differences between the fit and 
the data for the two models (note the different vertical scales). 



D. Self- force input 

To leading order, our analysis requires knowledge of 
L = F v along circular orbits with 3 < r < 6. In practice, 
we calculate E — —F t , which on circular orbits is related 
to L by L — E/Cl. To obtain E we used the self- force 
code developed in Ref. flfjj , which is a frequency-domain 
variant of the original Barack-Sago code @. The code 
was originally developed to deal with stable orbits, but it 
can handle unstable (circular) orbits without any further 
development. The code takes as input the orbital radius 
r and returns the value of E along that orbit. We have 
obtained accurate data for a dense sample of r values in 
the range 3 < r < 6. The numerical data and details of 
their derivation are given in Appendix [21 

Our numerical results can be approximated by the 
semi-heuristic fitting formula 

E ~ -rjr~ 5 z~ a (oo + a-yz + a 2 z 2 + a 3 z 3 ), (61) 

where z := 1 — 3/r and the parameters 

a ~ 1.773163, a ~ 2.87683, ai ~ -4.01414, 

a 2 ~ 10.5371, a 3 ~ -3.79087 (62) 

have been determined by a least-squares fit of E versus 
z. The r -5 decay at large r is well-known. We have 
modelled the blow-up of the self-force at the light ring 
by a single power z~ a , although we have no rigorous 
theoretical argument for this. The numerical data we 
have fitted to are also not very close to the light ring, 
with the smallest value of r at 3.02, corresponding to 
£7 ~ 4.15. 

A tentative theoretical argument for determining a is 
that at sufficiently high energy the metric perturbation 



should be proportional to the total energy E of the par- 
ticle and therefore the fluxes as seen by an observer at 
rest should be proportional to E 2 , so that 

E - iE 2 - z~ 3/2 . (63) 

[Note that, from dTOj) and ([27]). i = z~ 1/2 on circular 
geodesies.] A least-squares fit to ([oTj) where a = 1.5 
is held fixed results in a maximum fitting error of 17%. 
By comparison, the maximum fitting error with (|62D is 
0.5%, which is much smaller and comparable with the 
maximal estimated error of the self-force calculation (see 
Appendix \K§. The numerical self- force results and both 
fits are shown in Fig. [21 For our computations and plots 
in the remainder of the paper, we shall use the fit (|61[) 
with (j62p for the range r > 3.02 for which we have self- 
force data, corresponding to E < 4.15. 

Fig. S] shows the energy absorbed by the large black 
hole as a fraction of the total energy loss (radiated to 
infinity plus absorbed by the large black hole). This fig- 
ure is relevant for computing the total energy radiated 
to infinity, but not for our construction of the critical 
solution, where only the total energy loss matters. 



E. Integration of the critical solution 

In the presence of radiation reaction, for a given mass 
ratio rj, there is a single critical solution where E, L, f2, r, 
r and if are all monotonous functions of one another, with 
E and L decreasing from oo to finite values at plunge, r 
increasing from 3 to 6, and f2 decreasing from 3~ 3 / 2 to 
6 -3 / 2 , all as r and tp increase. In this paper we construct 
the critical solution only to leading order in 77. Then E, L 
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=l-3/r 



FIG. 4. The fraction of the total energy loss E/n that is 
absorbed by the large black hole, as a function of z := 1 — 3/r. 
This is based on the data in Appendix|X] Extrapolation gives 
a fraction of ~ 42% at the light ring itself. 



and f2 and r are known functions of one another, related 
by (|27M29[) . Furthermore, the ^-dependence of the critical 
solution is of the simple form r = r(t) = r(tp), where we 
define the "slow angle" 



ip := rjtp 



(64) 



in analogy to the slow time f . Finally, the ODEs obeyed 
by f(r) and <p(r) are separable, and can be solved by 
integration. 

We obtain tp(r) by (numerical) integration of 



dtp 
dr 



dip dt dr dE 
dt dr dE dr 



nt 



dE_ij_ 

dr E 



l r -B/a (r _3)-»(r-6)|, 



(65) 



after substituting (|5T|) . A very good closed form approx- 
imation to this integral in terms of elementary functions 
can be obtained by approximating (dtp / dr) / (r — 6) as a 
series in (r — 3) with three terms, multiplying back by 
(r — 6), and integrating this. The full result, the closed 
form approximation, and the ultrarelativistic approxima- 
tion discussed below are shown in Fig. [5] For fixed 77, this 
graph gives us the shape of the critical orbit. 

We only have reliable self-force data for r > 3.02, cor- 
responding to E < 4.15, but we are interested in the 
critical solution up to E = 00. We have therefore extrap- 
olated our fit (|6"1"T) with (|62[) down to r = 3 in producing 
Fig. [5j and have used this extrapolation to set <p(i) = 
at E = 00 as a matter of convention. Note, however, 
that almost the entire range of tp in this plot occurs for 
r > 3.02, where we do have reliable self-force data, with 
tp(3.02) ~ 0.07 compared to <p(6) ~ 2.6. This qualita- 
tive feature is robust, requiring only that E diverges at 
the light ring faster than (r — 3) _1 , that is a > 1. Our 
self-force data for r < 3.02 support this qualitative fea- 
ture robustly, with a ~ 1.77, while our crude theoretical 
model suggests a — 3/2. 




FIG. 5. Slow angle tp as a function of radius r in the critical 
solution. We have fixed tp = at r = 3. The full line repre- 
sents a numerical integration of (|65[) , with E(r) given by the 
fitting formula (|61|l . The dashed line represents the closed 
form approximation discussed in the text below Eq. (I65|l . 
and the dotted line the ultrarelativistic approximation ((69}. 
Note that here and in the following figure, the narrow section 
r < 3 < 3.02 of the domain is derived from an extrapolation 
of our self-force data. 



Hence, the total number of orbits in the critical solu- 
tion, from infinite energy down to plunge, is given by 



That 



<p(6) _ 0.41 

27T?7 



V 



(66) 



and we expect the numerical prcfactor to depend only 
weakly on our extrapolation of the self-force data. 
We obtain f(r) by (numerical) integration of 



dt 
dr 



j] dE 
E dr ' 



(67) 



A reasonably good closed form approximation for f(r) 
can be obtained by expanding (dt / dr) / (r—6) into a series 
in (r — 3) with four terms, multiplying back by (r — 6), 
and integrating. The full result and its approximations 
are shown in Fig. [6] Our self-force data suggest that f, 
too, is finite at infinite energy, and we use the convention 
that f = at E = 00. Again, compare f (3.02) ~ 0.02 
with f (6) ~ 9.62. 



F. The ultrarelativistic approximation 

In order to make quantitative statements about the 
critical solution beyond the energies covered by our self- 
force data, we shall assume that as E — > 00, E diverges 
like the leading order of (|6"Tj) . that is 



E 
V 



)Q — 5 



a (r - 3)~ 



(68) 



We shall refer to calculations where we only use the lead- 
ing power of (r — 3) throughout as the ultrarelativistic 
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FIG. 6. Slow time f as a function of radius r in the criti- 
cal solution. We have fixed f = at r = 3. The full line 
represents a numerical integration of (|67[) . The dashed line 
represents the closed form approximation for f (r) discussed 
in the text below Eq. (|67p and the dotted line the ultrarela- 
tivistic approximation (|70l) . 



approximation, and in this context we shall treat ciq > 
and a > 1 as undetermined. 

Combining the leading order of (|1)5|) with (p5)) and in- 
tegrating, we obtain 



37/2 



2a (a- 1) 



-3) c 



(69) 



Similarly, taking the leading power of (r — 3) in (|57|) . we 
obtain 



f(r) 



39/2- 



a (2a-l) 



(r-3r 



(70) 



In the ultrarelativistic approximation we can combine 
((30)) with (|69]) and (JTOj) to express r - 3, f, and i? 
all as powers of each other. In particular 



where 
C E := 



E(f) ~ C B f-^, 
- 5 a (2«-l)]- /3 , 



/3 



1 



2a 



which will be used later. With 

in the ultrarelativistic approximation, we also have 



C E 



-1-/3 



c 



1//3 



V3(l-/3) ~V3(l-/3) 
In the ultrarelativistic limit (|40[) becomes 

L~ E ~ C E ' 
where we have used (|71[) in the last equality. 



(71) 
(72) 

(73) 
(74) 
(75) 



G. Validity of the self-force and adiabatic 
approximations 

The self-force approximation to the equations of mo- 
tion is valid when the total energy mE of the small ob- 
ject is much smaller than the total energy M of the large 
black hole, that is, for 



E < if 



(76) 



Because the self-force diverges as E — > oo at fixed r/, we 
must also estimate where the adiabatic approximation 
breaks down at high energy. Hcuristically, we assume 
that the adiabatic approximation is valid if E changes at 
most by some small fraction i5 < 1 per orbit, or 



dlnE 5 
dip ' 2irr] 



(77) 



(The explicit rj on the right comes from the hat on the 
left.) From ([74]) . we find that in the ultrarelativistic 
regime this is equivalent to 



(78) 



We see that the adiabaticity condition ([75)) always im- 
plies the condition (fT6"|) for the self-force approximation 
itself to be valid if a > 3/2. (Recall our theoretically 
motivated value of a is 3/2, and our fitted value is 1.77.) 

Note that the adiabatic approximation breaks down 
also as the plunge is approached, because the effective 
potential becomes shallower, leading to a rapid increase 
in n at finite F r in (EFJ fnl. 



V. CRITICAL SCALING IN THE SELF-FORCE 
APPROXIMATION 

A. Perturbations of the critical solution 

We now consider a small perturbation of the critical 
solution that results from imperfect fine-tuning of the 
initial data, as discussed above in Sec. |TT1 namely 

r = r*(f) + <5r(r), L = i*(f) + SL(t), (79) 



where r* and L» are given by the adiabatic ansatz (|48I49[) 
above, while Sr and SL are just functions of r, that is, 
they are assumed to remain "fast" even in the limit of 
weak radiation reaction. The self-acceleration is per- 
turbed as 



F a = F2 + 0(rj)0(6r,6L), 



(80) 



where the factor 77 arises simply because any self- 
acceleration scales as 0(rf). Substituting this into the 
equations of motion and subtracting the critical solution, 
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we formally obtain 
1 



[V rr -{L , r )Sr + V, rL (L , r Q )SL] 



+0(Sr 2 ,6L 2 ,5rSL) 
SL = 0(r))0(5r, SL). 



0( V )0(6r,6L), 



(81) 
(82) 



We now neglect the error terms above, the first set as 
they are nonlinear in (Sr, SL), the second set as they are 
small corrections to the coefficients of a linear differential 
equation for (Sr, SL). 

In identifying the perturbed solution with the back- 
ground critical solution, we can adjust the origin of r in 
the perturbed solution. (This is a remnant of the usual 
gauge-dependence of spacetime perturbations.) As L»(f ) 
is not constant, wc can use this to set SL = at any one 
time, and hence, by SL ~ 0, at all times. Note that the 
neutral mode SL = constant present in the geodesic ap- 
proximation, which linked neigbouring critical solutions, 
has become a gauge mode in the presence of radiation re- 
action, where there is a single critical solution, and that 
we have used a simple form of gauge fixing to eliminate 
this mode. From now on we denote 8r in this gauge by 
x. We are left with the single equation of motion 



r(f)- 2 x, 



(83) 



which is equivalent to (|3T|) for the geodesic case, but 
where now 



T(f) :=T[L (f)} 



(84) 



is time-dependent. 

Rewriting the equation entirely in terms of f , we have 



(85) 



where as before a prime denotes d/dr. As r\ is small, we 
can use a WKB approximation. In terms of the "WKB 
time" 



T(r) 



dr 

r(f)' 



(86) 



where for definiteness we let T(0) = 0, the first-order 
WKB approximation to the growing and decaying solu- 
tions is 



•'■'4: 



(87) 



This is a good approximation for |T| <C 1, when the pref- 
actor varies much more slowly than the exponential. The 
functions T(f) and T(r) are shown in Figs. [JJ and [5] 

In the ultrarelativistic approximation, from (|75p . we 
have 



C E 



V3(l - P) 



\/3(l-/3)" 



ip, 



where we have used ([74")) in the last two equalities. In the 
same approximation, 



f ~ V3/3C 1 E - 2a i 1 E 2a - 



(89) 




10 



FIG. 7. Log-log plot of T(f) [defined in (|40], based on (|61l62f) 
(solid line), and the ultrarelativistic approximation Q75p (dot- 
ted line). (The divergence of V comes from the shallowness of 
the potential at r = 6.) 




FIG. 8. Log-log plot of T(f) [defined in J86], based on (pl62)l 
(solid line), and the ultrarelativistic approximation ()88|1 (dot- 
ted line). 



Hence ITI < 1 for 



E < if 



(90) 



which has the same power of rj as the criterion (|78|) for 
the critical solution to evolve adiabatically. Clearly, with 
6 < 1, flZHJ) implies ([90]) and hence, for a > 3/2, also 
([76")) . Hence the only criterion we need to impose is ([75)1 
(the critical solution evolves adiabatically), and we can 
then always use the WKB approximation to the critical 
solution. 

For r small enough (E large enough) that |r| 3> 1, the 
WKB approximation does not hold. However, we are 
then necessarily in the ultrarelativistic limit such that 
the approximation ([75]) for r(f) holds, and with this, (|53")l 
can be solved in closed form in terms of Bessel functions. 
As this result is not directly relevant for our main result, 
we give it in Appendix [Bl 
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B. Dependence on the initial data 

For definitencss, we parameterise the space of initial 
data (modulo a shift in time) by the initial energy at 
infinity Ei and the impact parameter b. We consider a 
generic near-critical solution in the whirl phase r, < r < 
Tf where it is close to the critical solution, so that 



x(t) ^A{E u b) 



B{E h b) 



x+{n)' 



(91) 



The factors x±(ti) have been introduced for later conve- 
nience. 

A and B are smooth functions of Ei and b, and by 
definition B vanishes for critical initial data. Therefore, 

A = A (Ei) + O(Sb), B = B 1 (E l )5b + 0(Sb) 2 , (92) 

where 



5b:=b-h(Ei,r)). 



(93) 



[The critical value b*(Ei,r]) of the impact parameter is 
known in closed form only for 77 = 0, when it is given by 
(|3"3"|). but its value does not matter in the following.] Wc 
must have A > and B\ > because x is positive and 
decreasing during the approach to the critical solution, 
and then becomes negative in the capture case b < &* , or 
returns to positive values in the scattering case b > 

Let the beginning and end of the whirl phase be defined 
by x(ji) = Xi and x(jf) = Xf, where Xi and Xf will be 
characterised later. If the whirl phase is sufficiently long, 
we can neglect the growing mode at t% and the decaying 
mode at Tf, so that 



A s Xi, 



Bn 



X 



x+jn) 

X+{Tf)' 



(94) 
(95) 



Note that Xi is always positive, while Xf is positive for 
initial data that scatter and negative for initial data that 
plunge. 

In the general discussion of critical scaling in Sec. [H] 
we introduced two constants C\ and C2 that were treated 
as unknown, and combined into a constant C3 in the final 
result ([4]). The equivalents of C\ and C2 in the current 
specific problem are B\ and Xf, respectively, and we shall 
now see that under certain conditions we can determine 
both, so that there is no undetermined constant left in 
our final results. 

In order to determine B\ we can use an energy argu- 
ment if we approximate E as constant during the ap- 
proach phase. This is justified when the energy loss dur- 
ing the approach phase is small compared to the energy 
loss during an extended whirl phase. We approximate 
(|87|) over a short time interval (covering the transition 
from approach to whirl) as 



x±{t) 



,±(T-rO/r(T„) 



x±(n). 



(96) 



Perturbing (fl"2"|) . which holds for geodesies, with respect 
to b, we obtain 

E 2 (L, b) - E 2 (L) ~ x 2 - T- 2 x 2 ~ -AT- 2 AB, (97) 



where the last equality follows from (|96l) . This equation 
becomes exact in the geodesic limit, where E is conserved 
and (|9"6")1 is exact. On the other hand, when E is con- 
served, from (|25[) we also have 



2L 



5b 



E 2 (L,b) - Ei(L) ~ —gf-Sb ~ -2[Et(L) - 1]-. (98) 

Note the factor E 2 - 1 in dE/db. It comes about because 
geodesic orbits with E = 1 are at rest at infinity, and 
so the impact parameter is not defined. More concretely, 
displacing a particle at rest at infinity in a direction or- 
thogonally to the line of sight to the large black hole be- 
fore dropping it corresponds to a rotation of the binary 
system and makes no difference to the infall. Because 
of this factor, some of our following results will be for- 
mally singular at E = 1. This is only because b becomes 
a bad parameter of the initial data. In practice, we are 
interested in initial data with E > 1 (and in particular 
B> 1), when this is not a problem. 

We now approximate E*(L) = Ei in (f9"5|) and T ~ 
r[L*(i£j)] =: Ti in ([9"7| because we approximate E and 
L as constant for r < t, and because Li ~ L^Ei) by 
fine-tuning. Hence we find 



AoBi 



r 2 (E 2 -i) 
2K 



(99) 



From (|94[) . (j95|) and ((99)) . and taking the logarithm for 
later convenience, we now obtain 



In 



\Sb\ 



In 



r?(£? - 1) 



2XiX 



In- 



-(r/) 



x+in) 



(100) 



This is our master equation for critical scaling: in princi- 
ple, it determines Ef in terms of the mass ratio rj and the 
initial data (Ei,b). In order to obtain a definite result, 
however, we still need to fix Xi and Xf. As we discuss 
below, there is a natural criterion for fixing them in the 
radiation reaction, but only a more arbitrary one in the 
geodesic limit. Moreover, the criterion we choose in the 
radiation reaction case has a singular limit as 7/ — > in 
the geodesic limit. Wc therefore discuss both cases sepa- 
rately, starting with the geodesic case. 



C. Critical exponents in the geodesic case 

The critical solution is by definition as circular as pos- 
sible in the presence of radiation reaction, sitting approx- 
imately at the top of the potential. Hence one possible 
criterion for Xi and Xf is that the radial velocity be 
given small fraction C e (for "eccentricity") of the tan- 
gential velocity, or 



x 



^ = C e . 



(101) 
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Note that C e plays the role of C2 in our general discus- 
sion in Sec. HH and while it must be small, its choice is 
arbitrary - we have it introduced it here because it is 
more intuitive than Xij. We have 



^E 2 -V \X h 



s\ 



X 



i,f\ 



r- l L Tr~ 1 L ^3TE' 



(102) 



where we have used (fTTj) and (fl~2|) in the first step, (|40l) in 
the second step, and the third approximation holds only 
in the ultrarelativistic limit. 

In the geodesic case, where the background solution is 
time independent, 

x± = e ±T ' v (103) 
holds exactly. Using this and (|100H102j) . we obtain 



Tf - r, 

were T is given by P0|) . Here 

E 2 -1 4 



-rin^i -TlnCe. (104) 
0* 



2C 2 r~ 2 L 2 2C 2 r 6C 2 



7^5 > (105) 



where the last approximate equality holds only in the 
ultrarelativistic limit r ~ 3. Clearly, C e inherits the ar- 
bitrariness of C e . The result (|104[) can be expressed in 
terms of orbits as 



1 I* 55 ' 1 r 
-7 In — 7lnC e 

0* 



where 



^ T 

7:= 2^ r 



2,1/2 



2ttL 



1/2 



(106) 
(107) 

(108) 
(109) 



In the ultrarelativistic limit, we have 

E 1 
1 

7 "27- 

The critical exponents T and 7 for the geodesic limit 
were previously obtained by Pretorius and Khurana 
for equatorial geodesies of Kerr. The additive constant 
parameterised by C e could only be determined more ex- 
actly by a full modelling of the perturbed orbit to replace 
the artifical split into a zoom-in, whirl and zoom-out or 
plunge phase. 

As an extension of the result for the time spent 
whirling, we can estimate the energy loss during the whirl 
phase as 



AE~(T f -n)E[r (L)}, 



(110) 



where (jf — Tj) is given in terms of the initial data by 
(|104jl . This approximation, where we neglect the radi- 
ation reaction on the orbit itself, is therefore consistent 
only for AE < £V 



D. Generalised critical scaling in the radiation 
reaction case 

With radiation reaction, a less arbitrary definition of 
Xi and Xf is given by equating the rate of change of 
the perturbation x with that of the background solution, 
that is 



(in) 



Furthermore, we have seen above that we can always use 
the WKB approximation (]87|) to the perturbation modes 
when the adiabatic self-force approximation to the back- 
ground critical solution is valid. Substituting our new 
definition (|111[) of Xi j into the left-hand side of the scal- 
ing master equation (|100[) , and the WKB approximation 
([57)1 into its right-hand side, (|100[) becomes 



, \Sb\ , T?T f (E?-l) Tf-Ti 1 T f , s 
In V- 1 -In \ t , ^ — + ~ In — • (H2 



1 0i' Of 

We can re-order this as 

_ In MM ^ If_Il + 2 In t? + K l + K f , (113) 
where we have defined the shorthands 

rf 2 (£f-i) 



K, := In ■ 



K f := In 



1r' 



^3/2 
/ 



(114) 
(115) 



Here we consider T, and Kf as functions of f (or 
any of the other slow variables E, r, ip). The only 77- 



dcpcndcnce is then the one explicitly shown in (| 1 13|) 
Alternatively, we can re-order f| 1 13[) as 



V 



K f 



In 



\Sb\ 



21nry. 



(116) 



This explicitly gives a known function of f/ (or Ef), on 
the left-hand side, as a function of the mass ratio 77 and 
the initial data (Ei,bi), on the right-hand side. This is 
the closest we can get to the fully explicit result Ef — 
Ef(rj,Ei,b) that we would ideally like. [We do not have 
it only for the technical reason that we cannot invert the 
function Tf/rj + Kf for Tf (or Ef) in closed form.] 

We can make (|116[) a bit more explicit, in terms of 
only powers and logarithms, when the entire whirl phase 
is ultrarelativistic, that is Ef 1. We have already 

given the ultrarelativistic approximations for T(f ), T(E) 
and T(ip) in (|88p . and we similarly find that that in the 
ultrarelativistic limit 



K ^ -2 In 2 



3/3 

In On 

2 



- +6/3 ) ln3 



(117) 



13 



K, 



Ki,Kf,ki T 




0.001 



FIG. 9. Log-linear plot of Ki(f) [defined in (|114| . based on 
(|61I62|I (solid line), and the ultrarelativistic approximation 
(|117[) (dotted line). The full curve ends because of the factor 
E 2 — 1 in the exact expression, which is approximated as E 2 
in the ultrarelativistic expression 




FIG. 10. Log-linear plot of Kf(f) [defined in 1)115] . based on 
(]61I62[I (solid line), and the ultrarelativistic approximation 
(PH (dotted line). 



-In 2 

P 



2(3 In 3 



In ao 



(118) 



[These can be readily re-expressed in terms of tp, E, or 
(r - 3) using ([oWTj) .] We have plotted f(f) above in 
Fig. |U and we plot K^f) and K f (r) in Figs. [51 and [TUI 
all based on (]61l62p and each compared against its ultra- 
relativistic approximation. Fig. [TT]plots all three against 
E instead of f . 

Note from Fig.QTJthat lnT, Ki and Kf vary compara- 
bly rapidly along the critical solution, but it is T /rj (with- 
out the logarithm), Ki and Kf that appear in our scaling 
result (|113|) . Hence one can suppress Ki and Kf in (|113j) 
as logarithmic corrections relative to the T terms. Doing 
this, suppressing also the constant term 2 In?;, and taking 




In E 



FIG. 11. Plot of (from top to bottom at large E) Ki, Kj and 
lnT (note the logarithm), all versus lni?, based on ()61l62p 
(solid lines), and the corresponding ultrarelativistic approxi- 
mations (dotted lines). 



the ultrarelativistic limit for clarity reduces (|113|) to 



-In- 



Tf - Ti 



C E 



V3(l - p) V / 



1-13 1-13 
Tf - Ti 



(119) 
(120) 



C 



1//3 



V3(l - P) V 5 
ff ~ ft- 



E l-w_ E }-ve\ 121) 



(122) 



It is interesting to note that the amount of fine-tuning of 
the initial data (measured by In \Sb\) required to achieve 
a given number of orbits on the critical solution is in- 
dependent (within the self-force approximation) of the 
mass-ratio rj. By comparison, better fine-tuning (smaller 
\Sb\) is required to achieve a given time on the critical 
solution, or a given amount of energy radiated, as r\ be- 
comes smaller. These scalings with r\ are not immediately 
intuitive. 

Finally, to see how the self-force result is related to 

The second 



the geodesic limit, we go back to (1112 
term on the left-hand side comes from the criterion (|111[) 
for determining Xi and Xf, and is simply not defined 
when the self-force is neglected. (Formally, it diverges 
as rj —> 0.) In the geodesic case we therefore adopted 
the more ad-hoc definition (|101|) of Xi and Xf. Hence, 
in the geodesic limit the left-hand side of (|112|) becomes 
— In |(56/6*| plus some constant related to Xi and Xf but 
independent of 5b. On the right-hand side of (]112)1 . we 



note r. 



r, 



r, so the second term vanishes, and 
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from (I55|) we have 



Hence we recover (|104[) as the 77 — ^ limit of (|112[) . but 
with the constant C e undetermined by (|112|) . 

VI. CONCLUSIONS 

Pretorius and Khurana Q suggested equatorial zoom- 
whirl orbits in Kerr as a toy model for the critical phe- 
nomena at the threshold of immediate merger in equal 
mass black hole binaries that they had observed in nu- 
merical relativity simulations. Here, we have adapted 
this toy model to make quantitative predictions in the 
limit where one of the two black holes is much larger than 
the other. We can then model the smaller black hole as a 
point particle of mass m in a background spacctime given 
by the larger one (assumed to be Schwarzschild with mass 
M = 1), and moving on an orbit determined by a grav- 
itational self-force proportional, to leading order, to the 
mass ratio r/ := m/M <C 1. 

In our model there is a single critical solution that 
evolves adiabatically along the sequence of unstable cir- 
cular geodesies, beginning with infinite energy and end- 
ing with a plunge when almost all the energy has been 
dissipated in gravitational waves. In this process, the 
critical orbit evolves from r = 3 and E — 00 (the light 
ring) to r = 6 and E = y/8/9 (the last stable circular or- 
bit), making only a finite number (approximately 0.41 / 77) 
of orbits. 

Any initial data fine-tuned to the critical impact pa- 
rameter are attracted to the critical solution, joining it at 
the appropriate energy. They leave it after a number of 
orbits determined by the degree of fine-tuning, to cither 
plunge or scatter. Beyond a certain degree of fine-tuning, 
they stay on the critical solution until it plunges. 

The type of critical solution we have described here is 
unfamiliar in general relativity, being neither self-similar 
(type II) nor stationary (type I). Rather, it avoids im- 
mediate plunge or scatter by remaining as circular as 
is compatible with radiation reaction. The appropriate 
symmetry ansatz is therefore a formal slow time expan- 
sion, in which the orbit evolves adiabatically from one 
circular geodesic orbit to an adjacent one under the ef- 
fect of energy and angular momentum loss. We might 
call this type la, for "adiabatic" . 

Our main result links the final energy (per rest mass) 
Ef of the small object to its initial energy Ei, the impact 
parameter b and the mass ratio r\. This is given in (|116l) . 
We have also given more explicit approximate expressions 
(neglecting subleading terms, and assuming for clarity 
that E f > 1) in Eqs. (|120til22[) . 

For the purpose of our calculation, there are four dif- 
ferent regimes for the initial energy: 

1) For \/8/9 < Ei < 4 we have accurate self- force data 
and so our results are reliable, but cannot be expressed 



in closed form. (We plot them, and the interested reader 
could reconstruct all our plots from the data and formu- 
las given in this paper.) Even in this regime, we had 
to extend previous self-force calculations much closer to 
the light ring at r = 3M, and it is clear that the self- 
force diverges there. It would be highly interesting in 
its own right to understand the origin of this divergence 
rigorously. 

2) For 4 < Ei < jy^ 1 /( 2Q - 2 ) j W e are extrapolating our 
numerical self- force data for E as (r — 3)~ Q , see Eq. (|68[) . 
(A best fit to our numerical data gives a ~ 1.77, but 
we also have a heuristic theoretical argument for a = 
3/2, which is a less good fit but not obviously wrong, 
see Fig. [3]) In this regime the motion is ultrarclativistic, 
which allows us to give explicit expressions. As 77 can be 
arbitrarily small, we give our (extrapolated) results up 
to Ei = 00. 

3) If a < 3/2, there is an the even higher energy regime 
? ^- 1 /( 2q -2) < E i < jy- 1 w hen the self-force approxima- 
tion is still valid, but the adiabatic approximation is not. 
(This regime is empty for a > 3/2.) It should be possi- 
ble to calculate in this regime in the foreseeable future, 
but at the moment of writing this, the gravitational self- 
force on non-geodesic orbits in the Schwarzschild spacc- 
time has not yet been calculated. (Similarly, it should 
become possible in the foreseable future to extend our 
self-force results to the case where the large black hole is 
Kerr.) 

4) Finally, for Ei > ?7 _1 , the (total, mostly kinetic) 
energy of the small particle becomes comparable to that 
of the large black hole, and the self-force approximation 
itself breaks down. (We come back to this regime just 
below.) 

Two surprising features of the critical solution at the 
threshold of immediate merger, at least in the extreme- 
mass ratio approximation, are its unfamiliar (approxi- 
mate) symmetry, and the fact that the fraction of E ra- 
diated per orbit increases with E fast enough that the 
critical solution evolves from E = 00 (or at least from 
E ~ rf 1 ^> 1, for any fixed finite mass ratio r\) to plunge 
at E ~ \/8/9 in a finite number of orbits. 

We conjecture that both these qualitative features hold 
also in the comparable to equal mass ratio case. Evidence 
that the fraction of energy radiated per orbit increases 
with energy is given by comparing the numerical relativ- 
ity results of [jf (1.0 - 1.5% at a boost of k ~ 0.22) and 
(~ 13% at a boost of k ~ 1.5 Q2). We should stress 
again that for any given mass ratio 77 (and spins) there 
is only one critical solution, which begins with infinite 
energy. We also conjecture that this solution radiates an 
infinite amount of energy over a finite, probably quite 
small, number of orbits. At low energies, an adiabatic 
approximation should be possible, where the critical so- 
lution evolves along a sequence of stationary solutions 
with a helical Killing vector under the influence of radia- 
tion reaction, but if our conjecture is correct this breaks 
down at sufficiently large energy. 

In the comparable mass case, the dependence of the 
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orbital radius of the critical solution on the energy will 
be very different from that in the extreme mass ratio 
case. For two point masses Mi := M and M2 := r\M in 
special relativity, where we set 77 < 1 by definition, with 
a relative boost k := —u\ui a > 1 towards each other, 
the total energy in their common centre of mass frame is 
given by 



(124) 



Furthermore, the hoop conjecture suggests that the crit- 
ical impact parameter for immediate merger is 



6* ~ EcM,i, 



(125) 



and this is borne out by numerical relativity simulations 
for comparable masses and large k [l3[. Finally, as the 
critical solution is by definition always between scattering 
and merger, we also expect 



(126) 



along the critical solution. In the limit krj ^> 1, (j 1 25[) 
gives us 6* ~ M \J2r\k and from (|126|) the critical solution 
spirals in from infinite orbital radius to the innermost 
stable circular orbit (ISCO). In the opposite limit kr\ <C 
1, (|125[) gives us 6* ~ M (a more precise result is of course 
b % ~ 3x^3 Al), the self- force approximation we have given 
here becomes valid, and the critical solution spirals out 
from the light ring to the ISCO. 

Constructing the critical solution presents a challenge 
at any mass-ratio: for self-force methods in the extreme 
mass-ratio case to reach higher energies and to go be- 
yond the adiabatic approximation to the orbit, and for 
effective one-body and numerical relativity methods in 
the comparable mass ratio case to construct the critical 
solution even at low energies. 
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Appendix A: Numerical data 

We present here the numerical data used to inform our 
"empirical" formula (|6ip for E(r). The data were gen- 
erated using the self-force code developed by Akcay in 
Ref. jiol ]. The code computes E in two different ways: 
(i) Directly, by calculating the component F t of the local 
gravitational self-force (per unit particle mass m) and us- 
ing E = —F t ; and indirectly, by numerically evaluating 
the flux of energy in gravitational waves radiated to infin- 
ity and down the event horizon, and using the energy bal- 
ance relation — E = (dEcw /dt)(dt/dr). Here dEcw/dt 



is the total flux of energy (per m) carried by the gravi- 
tational waves, and dt/dr = (1 — 3/r) -1 / 2 . dEQw/dt is 
made up of two pieces, dEQ W /dt and dEQ W /dt, respec- 
tively associated with radiation going out to infinity and 
absorbed by the black hole. The numerical evaluation of 
both F t and dE^^/dt are detailed in Ref. [l(j. We find 
that the two computations agree extremely well, differ- 
ing by no more than one part in ~ 10 10 for r > 3.2. At 
smaller radii the agreement is less good, but for none of 
the radii we considered (down to r = 3.02) is the dis- 
agreement greater than 1%. 

There is a practical limit on how close to r = 3 our 
numerical method can reach. Both above computations 
(of the local self-force and of the asymptotic fluxes) rely 
on a summation over multipole-mode contributions, with 



a practical cut-off at some I = l T 



Our calculations 



typically use Z max ~ 80 for the self-force and i max ~ 140 
for the fluxes. As the (ultra-relativistic) limit r — > 3 is 
approached, a beaming-like effect broadens the Z-mode 
power distribution and shifts it toward large I. For r — 3 
small enough it becomes computationally prohibitive to 
obtain a sufficient number of modes. In practice, we have 
not attempted to obtain accurate data below r ~ 3.02. 

Table U displays our numerical results for a sample of 
orbital radii in the range 3.02 < r < 6. We choose to 
show E data from flux calculations, which, we believe, arc 
slightly more accurate. For completeness we also show 
the flux dEQ W /dt absorbed by the black hole as a fraction 
of the total flux radiated. For circular orbits at the ISCO 
only very little (~ 0.3%) of the energy emitted goes down 
the hole. This fraction, however, becomes very significant 
for orbits closer to the light-ring, reaching as much as 
37.7% at r = 3.02. A rough extrapolation gives ~ 42% 
at the light ring itself. 



Appendix B: Perturbation modes of the critical 
solution in the ultrarelativistic limit 



In the independent variable T := rj T ', with f defined 
by (|86]), (83]) becomes 



d 2 x dlriT dx 
~dT i ~ dT dT 



x = 0. 



(Bl) 



The ultrarelativistic approximations (|75|) and (|88|) give 
T cx T /3 /( 1 -' 3 ), and so in the ultrarelativistic approxima- 
tion (IB1I) becomes 



d 2 x (3 I dx 
dr* ~ 1- (3TdT 

With the substitution 
x{T) =: T v y{T) cx f ^(T), v 



x = 0. 



2(1 -py 



(B2) 



(B3) 



we obtain 



T 2^V +T ^_ (T 2 + 2) (B4) 

dT 2 dT v iy ' v ; 
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r 


-rj^E 


% absorbed 


3.02 


83(1) 


37.7 


3.03 


39.9(2) 


36.3 


3.04 


23.62(5) 


35.1 


3.05 


15.64(1) 


34.1 


3.10 


4.232549(9) 


30.0 


3.15 


1.9223864(1) 


26.9 


3.20 


1.081977467(7) 


24.2 


3.30 


0.467163243(4) 


19.9 


3.40 


0.25012891219(4) 


16.3 


3.50 


0.1510176014(1) 


13.4 


3.60 


0.09864553444(1) 


11.0 


3.70 


0.06818392468(3) 


9.09 


3.80 


0.0491937318(2) 


7.50 


3.90 


0.03670796106(3) 


6.20 


4.00 


0.02814331203(5) 


5.15 


4.10 


0.02206146488(5) 


4.30 


4.20 


0.01761664969(8) 


3.61 


4.30 


0.01428856144(8) 


3.04 


4.40 


0.01174466751(8) 


2.57 


4.50 


0.00976536837(3) 


2.19 


4.60 


0.00820146283(7) 


1.88 


4.70 


0.00694902541(6) 


1.61 


4.80 


0.00593406093(6) 


1.39 


4.90 


0.00510284885(5) 


1.21 


5.00 


0.00441570494(4) 


1.05 


5.10 


0.00384285426(2) 


0.922 


5.20 


0.00336164181(3) 


0.810 


5.50 


0.00231155387(2) 


0.562 


5.75 


0.00173630600(1) 


0.424 


6.00 


0.00132984067(1) 


0.326 



TABLE I. Numerical data for the rate of energy loss, E — 
dE/dr, for a sample of sub-ISCO circular geodesic orbits of 
radius r. In the second column parenthetical figures show 
the estimated uncertainty in the last displayed decimals. The 
third column shows the percentage of energy radiated down 
the event horizon, 100 x (dE^ w /dt) /(dEaw /dt); all figures 
shown are significant. 



which is the modified Bessel equation with index v. 
Hence two linearly independent perturbation modes of 
the critical solution arc 



-(r) 



2tt 



7777(1-/?) 



#1/2 



V - 1 T{t) , (B6) 



with T(t) in these formulas now defined as the ultrarel- 
ativistic expression (|88[) . The normalisations have been 
chosen so that for large r such that |T| <C 1, these solu- 
tions have the asymptotic form ([571) . with T(f) defined 




'7(1 - P) 



(B5) 



0.001 



FIG. 12. Log-log plot of the growing perturbation mode x+(t) 
for n — 0.1 (lower curves) and 77 = 0.01 (upper curves). The 
solid lines are numerical solutions started up from (|B5|I at 
small r = 0.0001, the dashed lines are given by the WKB 
approximation (|87[) with the exact values of F(f ) and T(f), 
and the dotted lines are the ultrarelativistic (Bessel function) 
approximation (|B5|) . normalised to agree with the WKB ap- 
proximation. The power-law growth of both the exact mode 
and its ultrarelativistic approximation at small r visible here 
is in fact linear, as can be shown from (|B5[) . It can be seen 
that the WKB approximation is good at large f/small E, and 
the ultrelativistic approximation in the other regime. Note 
that our self-force data are valid only for f > 0.2, correspond- 
ing to r > 3.02. It is clear from the figure that starting the 
numerical integration with data from the ultrarelativistic ap- 
proximation at t = 0.2 instead of 0.001 would make little 
difference for n = 0.1 and very little difference for r; = 0.01. 



by and r(f ) defined by the ultrarelativistic expres- 
sion ([75)l . Hence we recover the WKB approximation 
(|87|) in the ultrarelativistic limit. On the other hand, as 
r — > 0, x + ~ f and .t_ ~ 1. Fig. IT121 and llUl show that 
an excellent approximation to the true growing mode is 
given by (|B5[) for small r and by (|87|) for large t, with 
the two approximations overlapping. 

In numerically evolving x+ towards increasing f from 
initial data given at small r by (|B5[) . we stably find the 
growing WKB solution at large r. Finding X- over its 
entire range is more difficult: in going forward from (|B6[) . 
numerical error triggers the growing mode, while going 
backwards from the WKB form (|87[) of X- results in a 
mixture of (|B5P and (|F36|) as f — > 0. Shooting from both 
sides would be required, but we have not done this. 
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FIG. 13. Log- log plot of the growing perturbation mode x+(t) 
for 77 = 0.01, showing the full range. The lines are as in 
Fig. 1121 From this and the preceding figure it is clear that the 
WKB approximation (I87|l is excellent for f > 0.2 for 77 = 0.1, 
and f > 0.001 for 77 = 0.01. Continuing this trend, the WKB 
approximation is good from smaller f (higher energies) for 
smaller 77, and it is always good essentially up to plunge. 
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